clear all;
clc;

% ----------------------- Geometry Parameters -----------------------------

% import data
patientgeom = load('femspine_patientgeom.txt');
cordgeom = load('patient_cordgeom.txt');

patnum = 5;

epidur  = 0; % 0 = intradural, 1 = epidural
subloc  = 1; % 1 = 1 mm above cord; 2 = 1 mm below dura
th_elec = 0; % 0, -10, -20

% location of electrode
if(epidur==1)
    eletag = 'epi_';
    fdist_elec = cordgeom(patnum,8);
else
    eletag = 'sub_';
    if(subloc==1)
        fdist_elec = cordgeom(patnum,9);
    elseif(subloc==2)
        fdist_elec = 1-cordgeom(patnum,9);
    else
        return;
    end
end
loctag = ['p',num2str(patnum),'_d',strrep(num2str(fdist_elec),'.','p'),...
          '_t',strrep(num2str(th_elec),'-','n')];

% cord geometry
delr_dura = 0.2;

betageom = patientgeom(patnum,:);
a_cord = betageom(1)/2;
b_cord = betageom(2)/2;
a_csf  = betageom(3)/2-delr_dura;
b_csf  = betageom(4)/2-delr_dura;
a_dura = a_csf + delr_dura;
b_dura = b_csf + delr_dura;
ycenter = 13.61;
csf_center = [0,ycenter];
ytopepi = 25;

y_cord = (ycenter-b_csf) + betageom(5)*2*b_csf;
cord_center = [0,y_cord];

% electrode location
if(epidur==1)
    ytmp = ycenter+b_csf+delr_dura;
    yelec = ytmp+(ytopepi-ytmp)*fdist_elec;
else
    yelec = (y_cord+b_cord)+(ycenter+b_csf-(y_cord+b_cord))*fdist_elec;
end
elec_center = [0,yelec];

th_rad = th_elec*(pi/180);
Arot  = [cos(th_rad),-sin(th_rad);sin(th_rad),cos(th_rad)];
rp = Arot*[0;yelec-y_cord];
xe = rp(1); ye = rp(2) + y_cord;
relec = [xe,ye];

% white matter
th_unit = 0:pi/20:2*pi;
rcord   = a_cord*b_cord./sqrt((b_cord*cos(th_unit)).^2+(a_cord*sin(th_unit)).^2); 
xcord   = cord_center(1) + rcord.*cos(th_unit);
ycord   = cord_center(2) + rcord.*sin(th_unit);
% csf space
rcsf   = a_csf*b_csf./sqrt((b_csf*cos(th_unit)).^2+(a_csf*sin(th_unit)).^2); 
xcsf   = csf_center(1) + rcsf.*cos(th_unit);
ycsf   = csf_center(2) + rcsf.*sin(th_unit);
% dura
rdura   = a_dura*b_dura./sqrt((b_dura*cos(th_unit)).^2+(a_dura*sin(th_unit)).^2); 
xdura   = csf_center(1) + rdura.*cos(th_unit);
ydura   = csf_center(2) + rdura.*sin(th_unit);

% grey matter
x1 = [0.00 ,0.62 ,1.24 ,1.87 ,2.09 ,2.24 ,1.45 ,2.03 ,2.22 ,2.03 ,1.21 ,0.62 ,0.00];
y1 = [-0.44,-0.12,0.96,2.04,2.06,1.88,0.43,-1.12,-1.62,-2.12,-2.18,-1.61,-0.88];
y1 = y1+cord_center(2);
xgrey = [x1,-fliplr(x1)];
ygrey = [y1,fliplr(y1)];
         
% ----------------------- Importing Data ----------------------------------

p_tag = num2str(patnum);

pol = -1;
noaxons = 200;

prefix = ['neg_',eletag,'bp_'];
suffix = [loctag,'.txt'];

cd(['patient',p_tag,'_thresholds']);
dcthresh = load([prefix,'dc_act_',suffix]);
[dcthresh(:,1),dc_si] = sort( dcthresh(:,1)*pol );
dc_xy = load(['p',p_tag,'_dc_xyloc.txt']);
dc_xy = dc_xy + ( cord_center'*ones(1,noaxons) )';
dc_xy = dc_xy(dc_si,:);

drthresh = load([prefix,'dr_act_',suffix]);
[drthresh(:,1),dr_si] = sort( drthresh(:,1)*pol );
dr_xy = load(['p',p_tag,'_dr_xyloc.txt']);
dr_xy = dr_xy + ( cord_center'*ones(1,noaxons) )';
dr_xy = dr_xy(dr_si,:);
cd ..;

cd(['patient',num2str(patnum),'_potentials']);
iunit = load([eletag,'bp_current_',loctag,'.txt']);
rtiss = 1/iunit/1e3;
cd ..;

% figure;
% hold on;
% for k = 1:200
%     plot(dr_xy(k,1),dr_xy(k,2),'k.');
%     gtmp = input('');
% end
% hold off;
% return

% ----------------------- Spatial Activation ------------------------------

xx = a_dura + 1.5;
a_bnd = [-xx,xx,ycenter-xx,ycenter+xx];

minth = min(dcthresh(:,1));
maxth = max(dcthresh(:,1));

frac_lv  = [0.35,0.4,0.45,0.5,0.55,0.6];
vstim_lv = minth + frac_lv*(maxth-minth);
nolev    = length(frac_lv);

figure;
for k = 5%:nolev
    
    %subplot(2,3,k);
    hold on;
%     plot(xcsf,ycsf,'k--','LineWidth',2);
    
    plot(relec(1),relec(2),'bo','MarkerSize',15,'MarkerFaceColor','b');
    imax_dc = round( frac_lv(k)*noaxons );
    ck_dr = ( drthresh(:,1) < dcthresh(imax_dc,1) );
    plot(dc_xy(1:imax_dc,1),dc_xy(1:imax_dc,2),'g.','MarkerSize',25);
    plot(dr_xy(ck_dr,1),dr_xy(ck_dr,2),'r.','MarkerSize',25);
    plot(xdura,ydura,'k-','LineWidth',5);
    plot(xcord,ycord,'k-','LineWidth',5);
    fill(xgrey,ygrey,[0.5,0.5,0.5],'LineWidth',3);
    
    istim = dcthresh(imax_dc,1)/rtiss;
    ttl = [num2str(frac_lv(k)*100),' % DC at I_{s} = ',num2str(istim,2),' mA'];
    title(ttl,'FontSize',25);
    xlabel('x (mm)','FontSize',25);
    ylabel('y (mm)','FontSIze',25);
    set(gca,'FontSize',28);
    %axis(a_bnd);
axis square;
%     axis square;
    hold off;
    
end
return

L1 = 'Electrode';
L2 = 'DCs activated'; 
L3 = 'DRs activated';
legend(L1,L2,L3,'Location','NE');
